1. Field of the Invention
This invention relates generally to the field of seismic data processing. Specifically, the invention is a method for the estimation and removal of artifact noise from seismic data.
2. Description of the Related Art
Geophysical prospecting is the activity of searching for deposits of valuable minerals or hydrocarbons, such as oil and gas. Seismic surveys are used in geophysical prospecting to determine the location of potential reservoirs of hydrocarbon deposits in subterranean rock formations. Seismic surveys measure the travel time of seismic waves from the starting locations at seismic sources through reflections off the interfaces between layers in rock formations to the ending locations at seismic receivers. Seismic processing combines the many travel times between source and receiver pairs and attempts to remove all undesired noise in the recorded seismic signal. The desired result is an image of the locations of the rock formations.
At present, seismic data processing techniques result in seismic images that are degraded by noise, or by seismic artifacts. A seismic artifact is any distortion in the seismic data that can impede the ability to accurately estimate reservoir properties of interest from seismic data. An example of a seismic artifact is the interference by a shallower object in the seismic survey of a deeper target. This type of interference, in general, includes the effect of the overburden.
FIG. 1 shows an elevation view illustrating a model of a seismic signal with this type of artifact. The interfering shallow object 101 could be any type of near surface geologic feature, such as shallow gas (shown), gas filled channel complexes, hydrate beds, or salt domes. The interfering shallow object 101 is located between a seismic survey line 102 on the surface and a deeper subsurface target 103, such as a potential hydrocarbon reservoir. The seismic survey line 102 could be located on a surface of the water in a marine seismic survey, or on a surface of the earth in a land seismic survey. The seismic survey line 102 contains the seismic sources 105 and the seismic receivers 106.
These artifacts may arise in several forms. Artifacts can be caused by acquisition parameters, processing parameters, and geologic factors. An acquisition artifact may arise from survey design and execution irregularities and is often spatially periodic or semi-periodic. A geologic artifact may arise from formations above a target horizon having geologic variation or an irregular shape and is often spatially non-periodic. A processing artifact may arise from imperfect processing algorithms and parameter choices, and may be either periodic or non-periodic.
Artifacts with the same spatial periodicity as the acquisition geometry are termed acquisition-related. The spatial periodicity of acquisition-related artifacts can be modeled using conventional modeling techniques. Once the periodicity is identified, that is, the wavenumbers have been determined, several conventional filtering techniques are available to remove (or, more generally, to manipulate) the energy with the same periodicity as the artifacts. However, determining the exact magnitude of the artifacts is beyond current modeling technology. Therefore, removal of the artifacts is often not possible with conventional techniques.
Other types of artifacts are spatially irregular. For example, surface obstacles encountered during data acquisition or geological features of the overburden positioned above a target interval may cause irregularities or shadowing effects in the data. In general, contaminated seismic data can lead to a distorted interpretation of formation properties, which can in turn lead to missed prospects, dry holes, or uneconomic development wells. Thus, it is desirable to be able to remove the effects of these types of artifacts from seismic data. No conventional technique is available for filtering random or pseudo-random artifacts.
The demand for precision in the measurement of formation properties from seismic data is intensifying. Seismic amplitude analysis has long been a key component of successful exploration and exploitation strategies. Increasingly subtle variations in the amplitudes of the seismic data are now being analyzed to achieve detailed areal delineation, estimate pay thickness and net/gross ratio, determine details of the depositional environment, and predict reservoir fluid content, lithology, and migration pathways. However, even within a given seismic data set, acquisition geometry and processing steps can produce significant artifacts and cause variability in seismic amplitudes that are unrelated to the formation properties of interest. The evaluation, quantification, and removal of artifacts are critical to realizing the full potential of quantitative seismic attribute analysis. By developing and applying techniques to remove the artifacts (or perturbations caused by conditions in the overburden) from the seismic data, the accuracy of predictions based on seismic attributes can be improved, and the associated risk in subsequent exploration activity can be lowered.
Currently, artifacts are filtered from seismic amplitude data using two broad approaches: empirical and deterministic. Empirical approaches examine the data in order to identify and manipulate the portion of the signal that appears non-geologic. Such approaches can significantly improve the interpretability of amplitude maps that are contaminated by artifacts such as water bottom multiples. Unfortunately, empirical approaches are susceptible to making merely cosmetic changes to amplitude maps without sufficiently improving the accuracy of the measurements of formation properties. Deterministic artifact mitigation approaches systematically derive corrections based on the causes of the artifacts. These approaches can lead to accurate corrections if all of the true causes of the artifacts are known and correctable. This, however, is rarely the case.
Three empirical methods are commonly used to remove periodic amplitude artifacts from post stack seismic data. These three methods are the moving window average, the discrete wavelet transform low-pass filter, and the fast Fourier transform filter. The latter two filter methods attempt to remove or correct all of the energy at affected periodicities (wavenumbers), but sometimes only deal with a portion of this energy. A comparison of the resolution of these three empirical methods is presented in the following paragraphs after a model of a seismic signal with a periodic amplitude artifact is defined.
The simplest form of a seismic signal with a periodic amplitude artifact is an additive artifact. In one dimension, a signal along a seismic line with an additive artifact can be defined in its most general form by the equation
f(x)=g(x)+p(x),xe2x80x83xe2x80x83(1) 
where
x is the Common Depth Point (CDP) location,
f(x) is the recorded seismic signal,
g(x) is the true geologic signal, and
p(x) is the noise signal.
FIG. 1 shows an example in which a seismic signal may have an amplitude artifact such as may be characterized by Equation (1). The seismic signal f(x), recorded at receiver 106, consists of the seismic representation of the geologic environment at the target 103 plus the imprint from local artifact 101.
The target or true geologic signal g(x) is the seismic signal that would be recorded if the artifacts created by the acquisition, processing, or geologic imprint were eliminated. It is desired to either solve for or estimate this geologic signal g(x). The perturbation noise signal p(x) represents the artifact noise created by the acquisition, processing, or geologic imprint, and may be either random, semi-periodic, or periodic. The amplitude of perturbation p(x) is unknown, although the overall character of its spatial variation or periodicity is known or may be estimated empirically. Thus p(x) can be written as a product of an unknown amplitude modulation factor and the known or estimated spatial periodicity factor. Doing this, Equation (1) can be rewritten as
f(x)=g(x)+xcex5(x)*n(x),xe2x80x83xe2x80x83(2) 
where
xcex5(x) is the amplitude modulation factor of the noise signal, and
n(x) is the spatial periodicity factor of the noise signal.
The signal xcex5(x) represents the unknown amplitude by which the spatial periodicity factor n(x) is added to the true geologic signal g(x). Although referred to as the spatial periodicity factor, it is understood that n(x) may be either periodic or random, and as indicated above, may be known or estimated empirically. If an estimate of xcex5(x) were available, then the estimated solution of Equation (2) would be given by
ĝ(x)=f(x)xe2x88x92{circumflex over (xcex5)}(x)*n(x).xe2x80x83xe2x80x83(3) 
Here the hat or caret symbol ({circumflex over ( )}) denotes estimates of the signals g(x) and xcex5(x), respectively.
FIGS. 2A-2D show a model of a seismic signal with a periodic amplitude artifact, in each case plotted as amplitude versus CDP number. FIG. 2A is a plot of the amplitude of the recorded seismic signal f(x) 201 versus common depth point (CDP) number. FIG. 2B shows the spatial periodicity signal n(x) 202 of an acquisition artifact and FIG. 2C shows the amplitude modulation signal xcex5(x) 203 of the artifact. FIG. 2D shows the target or true geologic signal g(x) 204.
The first commonly used method for removing periodic amplitude artifacts is the moving window average. An analysis of the resolution of this method is shown in FIGS. 3A-3C, using the model shown in FIGS. 2A-2D.
FIG. 3A shows a comparison in which signal f(x) 301 from FIG. 2A and the spatial periodicity signal n(x) 302 from FIG. 2B are plotted on the same axis. The horizontal scale (CDP number) in FIGS. 3A-3C has been expanded, as compared to FIGS. 2A-2D, for clarity. FIG. 3B shows a comparison of the estimated geologic signal ĝ(x) 303 with the geologic signal g(x) 304 from FIG. 2D. Here, the signal ĝ(x) 303 is estimated using a moving window average. FIG. 3C shows a comparison of the averaged noise signal n(x) 305 with the averaged noise product signal xcex5(x)*n(x) 306. Note that {circumflex over (xcex5)}(x)*n(x)=f(x)xe2x88x92ĝ(x), the difference between the original recorded seismic signal and the final estimated geologic signal. The fundamental periodicity of the artifacts determines the proper window size in the moving window average. FIG. 3D shows the loss of resolution, a typical shortcoming of moving window averaging techniques.
The poor fit between the estimated geologic signal ĝ(x) 303 and the actual geologic signal g(x) 304 in FIG. 3B is captured by R2 the square of the correlation coefficient between the two signals R2=0.659. A perfect fit of R=R2=1 implies that 100% of the variance in the signal is captured. R and R2 are in non-dimensional units, always varying between xe2x88x921 and 1 or 0 and 1, respectively.
The poor fit between ĝ(x) 303 and g(x) 304 in FIG. 3B is also captured by the Root Mean Square Error, RMSE=228.8, given in the same units as the signal (here, the units of seismic amplitude). The Root Mean Square Error represents the average variability of the difference between the target signal g(x) and the predicted (or filtered) signal ĝ(x). A smaller RMSE indicates better results, since it indicates smaller differences. Such statistical measures will be used throughout the comparisons.
The second commonly used method for removing periodic amplitude artifacts is the discrete wavelet transform low-pass filter. The discrete wavelet transform low-pass filter method attempts a straightforward removal of energy at the fundamental frequency of the artifact and that frequency""s harmonics. An analysis of the resolution of this method is shown in FIGS. 4A-4C, using the model shown in FIGS. 2A-2D.
FIG. 4A shows the energy distributions of the input signal f(x) 401 and the spatial periodicity signal n(x) 402 before application of the discrete wavelet transform low-pass filter, while FIG. 4B shows the energy distributions of f(x) 403 and n(x) 404 after application of the discrete wavelet transform low-pass filter. FIGS. 4A and 4B are shown as plots of energy versus periodicity levels. As will be understood to those skilled in the art, wavelet analysis signals are decomposed into scales or levels of variability. For example, a signal of 128 points in length might be decomposed into k=7 scales (27=128) or k+1 levels. These scales, or levels, are analogous to bands of frequency, or wavenumber, in traditional Fourier spectral analysis. The highest scales are associated with the smallest (shortest) trend features. The lowest scales capture the longest trend features. For example, for a signal of 128 points in length, level=8 corresponds to a periodicity of 4 points (twice the Nyquist frequency); level 7 corresponds to an 8-point periodicity, and so on. FIG. 4C shows a comparison of the estimated geologic signal ĝ(x) 405 with the actual geologic signal g(x) 406 from FIG. 2D, as a plot of amplitude versus CDP number. Here, the signal ĝ(x) 405 is estimated using a discrete wavelet transform low-pass filter.
As with the moving window average method described above, there is a significant loss of resolution with the discrete wavelet transform low-pass filter method. The fit between the estimated geologic signal ĝ(x) 405 and the actual geologic signal g(x) 406 in FIG. 4C, as measured by R2=0.766 and RMSE=187.7, is only marginally better than for the moving window average technique as shown in FIG. 3B. In this example all the energy at scales 4 through 7 was removed, as indicated in FIG. 4B.
The third commonly used method for removing periodic amplitude artifacts from seismic data is the fast Fourier transform (FFT) filter. The fast Fourier transform filter method attempts to remove all the energy in the periodicities or wavenumbers affected by the artifacts. An analysis of the resolution of this method is shown in FIGS. 5A-5C, using the model shown in FIGS. 2A-2D.
FIG. 5A shows the wave number distributions of the input signal f(x) 501 and the spatial periodicity signal n(x) 502 before the application of a traditional fast Fourier transform filter, while FIG. 5B shows the wave number distributions of f(x) 503 and n(x) 504 after the application of the fast Fourier transform filter. FIGS. 5A and 5B are shown as plots of logarithm of the amplitude of the signal versus wavenumber. FIG. 5C shows a comparison of the estimated geologic signal ĝ(x) 505 with the actual geologic signal g(x) 506 from FIG. 2D, as a plot of amplitude versus CDP number. Here, the signal ĝ(x) 505 is estimated using a fast Fourier transform filter.
The energy in the wavenumbers contaminated by the artifact has been reduced, and the loss of resolution is less than in the examples shown in FIGS. 3B and 4C. However, the fit between the estimated geologic signal ĝ(x) 505 and the actual geologic signal g(x) 506 in FIG. 5C, as measured by R2=0.668 and RMSE=231.9, is only comparable, if not worse.
For either of the discrete wavelet transform or fast Fourier transform techniques to work, it is essential to have an estimate of the scales or periodicities of the artifacts. If the artifacts are spatially random (e.g., like those produced by geologic variations in the overburden), then any technique based on the discrete wavelet transform or the Fourier transform will give subjective or non-repeatable results as it fails to preserve the main characteristics of the noise.
The examples presented above in FIGS. 3 through 5 are typical of the results obtained by applying the three exemplified methods to one-dimensional signals and provide the context in which the impact of the method of the present invention can be evaluated. The examples show that there is a need for a relatively inexpensive filtering method that is able to identify and then remove the noise due to artifacts from seismic data without introducing bias or aliasing and without loss of resolution. Additionally, this method needs to be applicable to both periodic and non-periodic (random) perturbations.
The invention is a method for estimating and removing artifact noise from a seismic signal. First, the seismic signal is represented as a combination of a geologic signal and a noise signal. Then the representation of the seismic signal is decomposed into a linear combination of orthonormal components. One or more components are identified in which the geologic signal and the noise signal are uncorrelated. The noise signal is expanded as a product of an unknown noise amplitude modulation signal and a known noise spatial periodicity signal at the identified components. The expansion is solved for an estimate of the noise amplitude modulation signal at the identified components. Finally, the noise signal is estimated by multiplying the estimated noise amplitude modulation signal by the noise spatial periodicity signal. Additionally, the estimated noise signal can be removed from the seismic signal to generate an estimate of the geologic signal.